\documentclass[11pt]{article}
%% \topmargin -1.5cm        % read Lamport p.163
%% \oddsidemargin -0.04cm   % read Lamport p.163
%% \evensidemargin -0.04cm  % same as oddsidemargin but for left-hand pages
%% \textwidth 16.59cm
%% \textheight 21.94cm 
%$ %\pagestyle{empty}       % Uncomment if don't want page numbers
\parskip 7.2pt           % sets spacing between paragraphs
%% %\renewcommand{\baselinestretch}{1.5} % Uncomment for 1.5 spacing between lines
\parindent 0pt		 % sets leading space for paragraphs

\usepackage{fullpage}
\usepackage{listings} % For source code
\usepackage[usenames,dvipsnames]{color} % For colors and names
\usepackage[pdftex]{graphicx}
\usepackage{subfigure}

\usepackage{fixltx2e}

\definecolor{mygrey}{gray}{.96} % Light Grey
\lstset{ 
        language=[ISO]C++,              
        tabsize=3,                                                     
        basicstyle=\tiny,               
        numbers=left,                   
        numberstyle=\tiny,              
        stepnumber=2,                   
        numbersep=5pt,                  
        backgroundcolor=\color{white}, 
        %showspaces=false,              
        %showstringspaces=false,        
        %showtabs=false,                
        frame=single,                    
        tabsize=3,                          
        captionpos=b,                   
        breaklines=true,                
        breakatwhitespace=false,        
        %escapeinside={\%*}{*)},        
        commentstyle=\color{BrickRed}   
}

\title{Northeastern University \\
  Department of Electrical and Computer Engineering \\
  - \\
  ECE5667 \\
  Final Project}
\date{\today}
\author{
  Instructor: Deniz Erdogmus \\
  TA: Hooman Nezamfar \\
  Lab Partner: Andrew Lai \\
  - \\
  Author's Name: Paul Ozog
}

\begin{document}

\begin{titlepage}
  \maketitle
  \thispagestyle{empty}
\end{titlepage}

\tableofcontents

\pagebreak

\section{Introduction}
For our final project, we apply knowledge from the previous labs to implement and analyze a variety of digital filters on the ADSP-BF535 Ez-Kit Lite.
For Part 1, we will remove an constant-frequency interference tone from an audio signal using an Infinite Impulse Response (IIR) filter with the help of MATLAB's \texttt{fdatool}.
Part 2 involves deriving a filter analytically to remove an echoed portion of a speech signal.  For each part, we will modify the included \texttt{Talkthrough} program to implement the filters on the ADSP-BF535.

\section{Results and Analysis}

\subsection{Part I}
\subsubsection{Exercise 1: Frequency Analysis Using MATLAB}
To find the spectrum of the corrupted audio in Exercise 1, we used the MATLAB script shown below.  One should note that the audio in the \texttt{.wav} file is sampled at 11025 Hz.
\begin{verbatim}
[x, fs] = wavread('betho2n.wav');
fiveSec = x(100000:100000+fs*5.0);

mag = abs(fftshift(fft(fiveSec)));
freq = linspace(-fs/2, fs/2, fs*5+1); %for abs freq
magNorm  = mag(floor(length(mag))/2:end);
freqNorm = linspace(0, fs/2, length(magNorm))/fs; %for normalized freq

figure; plot(freq, mag)
xlabel('Frequency (Hz)');
ylabel('FFT Magnitude');
figure; plot(freqNorm, magNorm)
xlabel('Normalized Frequency (\pi rad/sec)');
ylabel('FFT Magnitude');
\end{verbatim}

This resulted the plots of the FFT magnitude, shown in Figure \ref{ex1}.  Clearly, the interference tone has a frequency of 1300 Hz.  

\pagebreak

\begin{figure}[h]
  \begin{center}
    \subfigure[Normalized Frequency] {
      \includegraphics[scale=0.5]{ex1-1.pdf}
      \label{ex1:1}}
    \subfigure[Actual Frequency] {
      \includegraphics[scale=0.5]{ex1-2.pdf}
      \label{ex1:2}}
    \caption{Spectrum of Corrupted Signal}
    \label{ex1}
  \end{center}
\end{figure}

\subsubsection{Exercise 2: Filter Design}
To design the filter, we used \texttt{fdatool} to derive the IIR coefficients.  The settings are shown in Table \ref{table:ex1}.  Once designed, the filter coefficients were exported to the Workspace using the \texttt{Export} feature.

\begin{table}[h]
  \caption{Exercise 1 IIR Filter Specs}
  \centering
  \begin{tabular} {c c}
    \hline\hline
    Specification & Value \\ [0.5ex]
    \hline
    Filter Type & Butterworth (IIR) \\
    Order & 2 \\
    \begin{math}
      f_c
    \end{math} & 13.0 kHz \\
    \begin{math}
      \Delta f
    \end{math}
    & .30 kHz \\ [1ex]
    \hline
  \end{tabular}
  \label{table:ex1}
\end{table}

This resulted in the Z-Domain coefficient vectors 
\begin{math}
  {\bf B}
\end{math} 
 (numerator) and
\begin{math}
  {\bf A}
\end{math}
 (denominator) as follows:

\begin{verbatim}
B = [1, -1.47597299627921, 1.00000000000000];
A = [1, -1.45115994760889, 0.966377367698636];
B_Gain = 0.983188683849318;
A_Gain = 1;
\end{verbatim}

and the performance plots are shown in Figure \ref{ex2}.  Note that the poles are very close to the unit circle, yet the system is stable overall.  However, it is important to consider quantization error in filter coefficients and calculation of difference equations when implementing this on the ADSP-BF535.

\begin{figure}[h]
  \begin{center}
    \subfigure[Magnitude Response] {
      \includegraphics[scale=0.5]{ex2-1.pdf}
      \label{ex2:1}}
    \subfigure[Pole-Zero Plot] {
      \includegraphics[scale=0.5]{ex2-2.pdf}
      \label{ex2:2}}
    \caption{IIR Interference Removal Filter}
    \label{ex2}
  \end{center}
\end{figure}

\pagebreak

\subsubsection{Exercise 3: MATLAB Simulation}

Armed with the interference removal filter's coefficients we could simulate the performance in MATLAB using the \texttt{filter} command like so:

\begin{verbatim}
y = filter(B*B_Gain,A*A_Gain,fiveSec);
\end{verbatim}

which resulted in the plots in Figure \ref{ex3}.  One should note that there is no spike at 1.3 kHz in Figure \ref{ex3:2}.  Therefore, we have removed the interference tone using MATLAB as a filter simulation environment. 

\begin{figure}[h]
  \begin{center}
    \subfigure[Normalized Frequency] {
      \includegraphics[scale=0.5]{ex3-1.pdf}
      \label{ex3:1}}
    \subfigure[Actual Frequency] {
      \includegraphics[scale=0.5]{ex3-2.pdf}
      \label{ex3:2}}
    \caption{Spectrum of Filtered Signal}
    \label{ex3}
  \end{center}
\end{figure}

\subsubsection{Exercise 4: Filter Implementation on DSP Board}
To implement on the ADSP-BF535, we had to redesign the filter to have a sampling rate of 48 kHz.  Every other specification remained the same.  This resulted in the following filter coefficient vectors 
\begin{math}
  {\bf B}
\end{math}
and
\begin{math}
  {\bf A}
\end{math}
with elements 
\begin{math}
  b_k
\end{math}
and
\begin{math}
  a_k
\end{math}
:

\begin{verbatim}
B = [1.0000, -1.97112731672697, 1.0000];
A = [1.0000, -1.96341695673344, 0.992176700177507];
\end{verbatim}

To implement a IIR filter on the DSP board, we had to convert the transfer function to a time-domain difference equation.  Since the Z-Domain transfer function was of the form:

\begin{equation}
  H(z) = \frac{Y(z)}{X(z)} = \frac{\sum_{k=0}^{k=N_z} b_k z^{-k}}{\sum_{k=0}^{k=N_p} a_k z^{-k}}
\end{equation}

We simply took the inverse Z transform.  Using algebra to express the current value of the filter in terms of previous values of the filter and input signal becomes:
\begin{equation}
  y[n] = \frac{-\sum_{k=1}^{N_z} a_k y[n-k] + \sum_{k=0}^{N_p} b_k x[n-k]}{a_0}
\end{equation}

These equation and coefficients were transfered to the \texttt{Process\_Audio\_Data()} function as follows:

\begin{verbatim}
void Process_Audio_Data (void) {
...
  CircBuffIn[i1] = sLeft_Channel_In;

  CircBuffOut[i1] = (b0*CircBuffIn[i1] + b1*CircBuffIn[i2] + b2*CircBuffIn[i3]
                     - a1*CircBuffOut[i2] - a2*CircBuffOut[i3]);

  sLeft_Channel_Out = CircBuffOut[i1];
  sRight_Channel_Out = CircBuffOut[i1];
...
}
\end{verbatim}

The the signal is mono, so only one circular buffer on the left channel is required to implement the filter.  Also, the array indexes are ensured to not be out-of-bounds. 

Before running the program, we ensured that the board's audio codec's input source selected \texttt{Line-In} and the jumper was set to the line-in mode.

Initially, the filter was unstable (the output quickly blew up to a loud meaningless signal).  We determined two causes.  First, the filter coefficients were overly quantized, so we simply added more decimal places to the coefficients.  Second, we initially computed the filter output using \texttt{short}'s.  Switching to floating point precision greatly helped to alleviate this problem.  A combination of these two improvements made the filter completely stable for all inputs.  

Connecting the function generator to the board input, and using a scope to analyze the output, we observed near-infinite attenuation at the center frequency of 1.3 kHz.  The scope display is shown in Figure \ref{fig:ex4}.

\begin{figure}[h]
  \begin{center}
    \includegraphics[scale=0.4]{PRINT_02.png}
    \caption{Time domain response when {f}\textsubscript{in} = 1.32kHz}
    \label{fig:ex4}
  \end{center}
\end{figure}

The output of the filter in Figure \ref{fig:ex4} is simply noise.  At all other frequencies outside the stopband, we observed insignificant attenuation of the input sinusoid.

\subsubsection{Exercise 5: Spectrum Analysis}

To examine the signal spectrum before filtering, we simply connected the computer's audio output to the input of a oscilloscope.  We provide two plots, shown in Figure \ref{fig:ex5:1}.


\begin{figure}[h]
  \begin{center}
    \subfigure[0 to 11025 Hz] {
      \includegraphics[scale=0.4]{PRINT_03.png}}
    \subfigure[300 to 2300 Hz] {
      \includegraphics[scale=0.4]{PRINT_04.png}}
    \caption{Spectrum of Unfiltered Signal}
    \label{fig:ex5:1}
  \end{center}
\end{figure}

Notice that there is a strong spike at exactly 1.3 kHz.  Its peak value is about 50.6 dB higher than the rest of the signal's frequencies.

Next, we fed the audio to the ADSP-BF535, and connected the output of our filter to the scope.  The spectrum of the result is shown in Figure \ref{fig:ex5:2}.  Clearly, the interference tone is removed from the signal because there is no spike at 1.3 kHz.

\begin{figure}[h]
  \begin{center}
    \subfigure[0 to 11025 Hz] {
      \includegraphics[scale=0.4]{PRINT_05.png}}
    \subfigure[300 to 2300 Hz] {
      \includegraphics[scale=0.4]{PRINT_06.png}}
    \caption{Spectrum of Filtered Signal}
    \label{fig:ex5:2}
  \end{center}
\end{figure}

\subsection{Part II}

\subsubsection{Exercise 1: Signal Representation}
A simple echo system with delay of 
\begin{math}
  \Delta t
\end{math}
seconds and a scaling of 
\begin{math}
  \alpha
\end{math}
is characterized by the following difference equation:

\begin{equation}
  y[n] = x[n] + \alpha x[n - \Delta t F_s] = x[n] + \alpha x[n - N]
\end{equation}

where x[n] is the input, y[n] is the output, and 
\begin{math}
  F_s
\end{math}
is the sampling rate in samples per second.  N is the closest integer to 
\begin{math}
  \Delta t F_s.
\end{math}


In the Z-domain, the system can be described by 

\begin{equation}
  H(z) = \frac{Y(z)}{X(z)} = \frac{1 + \alpha z^{-N}}{1}
\end{equation}

To determine the impulse response h[n], we used the following coefficient vectors in a MATLAB script:
\begin{verbatim}
A = 1;
B = zeros(1, N+1);
B(1) = 1; B(end) = alpha;
\end{verbatim}
These coefficients can be used for the standard selection of MATLAB commands like \texttt{impz} and \texttt{filter}.  The impulse response of h[n] is shown in Figure \ref{fig:2:ex1}.

\begin{figure}[h]
  \begin{center}
    \includegraphics[scale=0.6]{ex1_part2.pdf}
    \caption{Impulse response of h[n]}
    \label{fig:2:ex1}
  \end{center}
\end{figure}

\subsubsection{Exercise 2: Filter Design}
\label{foobar}
To remove the echo, we may use the following system with input 
\begin{math}
  Y(z)
\end{math}
and output
\begin{math}
  W(z)
\end{math}

\begin{equation}
  H_{inv}(z) = \frac{W(z)}{Y(z)} = \frac{1}{1 + \alpha z^{-N}}
\end{equation}

In the time domain, the output can expressed recursively as:

\begin{equation}
  w[n] = y[n] - \alpha w[n - N]
  \label{eq:2:diff}
\end{equation}
The impulse response of 
\begin{math}
  H_{inv}(z)
\end{math}
is shown in Figure \ref{fig:2:inv}.  It was plotted using the \texttt{impz} command with the coefficient vectors of the 
\begin{math}
  H_{inv}(z)
\end{math}
system:
\begin{verbatim}
impz(A,B);
\end{verbatim}

\begin{figure}[h]
  \begin{center}
    \subfigure[Inverse of h]{
      \includegraphics[scale=0.4]{echo-removal-impulse-paul.pdf}
      \label{fig:2:inv}}
    \subfigure[Cascaded System] {
      \includegraphics[scale=0.4]{ex2_part2.pdf}
      \label{fig:2:ex2}}      
    \caption{Impulse Responses of Echo Removal Systems}
  \end{center}
\end{figure}

The plot of the cascaded systems
\begin{math}
  H(z) H_{inv}(z)
\end{math}
is shown in Figure \ref{fig:2:ex2}.  

We expect a Dirac delta impulse response because clearly
\begin{math}
  H(z) H_{inv}(z) = 1.
\end{math}
The cacased system was determined using nested \texttt{filter} commands:
\begin{verbatim}
delta = zeros(1,100);
delta(length(delta/2)) = 1;
cascade_out = filter(A, B, filter(B,A,delta));
\end{verbatim}

This analysis using MATLAB shows that when 
\begin{math}
  x[n]
\end{math}
(the uncorrupted signal) is recoverable when passed through the cascaded system:
\begin{equation}
  \left( x[n] * h[n] \right) * h_{inv}[n] = x[n]
\end{equation}

For our group's \texttt{.wav} file, 
\begin{math}
  \alpha = 0.6,
\end{math}
and
\begin{math}
  N = 2000
\end{math}
(even though we had to use N = 2001 for the MATLAB simulation using \texttt{sound} to cancel the echo):

\begin{verbatim}
N = 2001; Fs = 7040; alpha = 0.6;
echoedSpeech = wavread('testing_2000_0.6.wav');
sound(filter(A, B, echoedSpeech), Fs)
\end{verbatim}

This resulted in producing a clear speech signal without any echo at all.  Clearly, the echo removal filter worked in simulation.

\subsubsection{Exercise 3: Filter Implementation}

To fit the much higher sampling rate of 48 kHz, the parameter N had to be adjusted:
\begin{equation}
  N = \left( \frac{2000}{7040} \right) F_s = 13643
\end{equation}
however the scale parameter
\begin{math}
  \alpha
\end{math}
did not change. 

The main lines of \texttt{Process\_Audio\_Data()} were:
\begin{verbatim}
...
i1 = ind;
i2 = ind - N;
i2 = (i2 < 0) ? i2 + ARRAY_SIZE : i2;

CircBuffIn[i1] = sLeft_Channel_In;
CircBuffOut[i1] = CircBuffIn[i1] - alpha*CircBuffOut[i2];

sLeft_Channel_Out = CircBuffOut[i1];
sRight_Channel_Out = CircBuffOut[i1];  
...
\end{verbatim}

Though the above code snippet looks convoluted, all it is doing is just implementing Equation \ref{eq:2:diff}. The result is stored in a circular buffer of 14000 samples, and is assigned to the output of the audio codec. 

To test the filter, we recorded a 20-second speech .wav file.  Using MATLAB, we applied the same length of echo that was given to our group's sample \texttt{.wav} file.  We connected the output of the computer's audio jack to the input of the board.  Speakers were connected to the board output, and we measured the results qualitatively. 

The echo removal filter worked, however when being filtered on the board the echo removal was not as clean compared to the MATLAB simulation.  When implemented on the ADSP-BF535, the echo wasn't completely removed.  If one were to listen carefully, very faint echos may still be heard through the speakers.  However, for all practical purposes the IIR filter still removed a sufficient amount of echo. 

\subsubsection{Sensitivity to Delay}
The filter When N (the delay in samples) was increased beyond
\begin{math}
  \pm 7
\end{math}
samples, the performance of the echo removal filter became unacceptable for a sampling rate of 48000 Hz.  This region corresponds to approximately 1 sample when the sampling rate is only 7040 Hz.  This explains why the MATLAB simulation in Section \ref{foobar}  worked for N = 2001, but didn't work at all for N = 2000.

\subsubsection{Sensitivity to Scale}
When the scale factor 
\begin{math}
  \alpha
\end{math}
was increased or decreased by .1, the effect on the quality of the echo removal filter was not noticeable.  If 
\begin{math}
  \alpha
\end{math}
were decreased from .6 to say .45, the filter would not remove a sufficient amount of echo.  On the other hand, if 
\begin{math}
  \alpha
\end{math}
were increased to .75, the filter would actually {\it produce} echo because the scale factor is so great that the signal's echo gets over-compensated.

\section{Conclusion}
Using a second-order Butterworth IIR, we were successfully able to remove the interference tone of 1.3 kHz in Part I.  From a qualitative perspective, the interference tone was completely gone and the filter did not affect the overall quality of the music.  For Part II, the implementation of the echo removal system sufficiently removed the echos from a corrupted speech signal.  We showed that to operate the properly, the filter must have a delay
\begin{math}
  N = 13643 \pm 7
\end{math}
samples, and the scale factor 
\begin{math}
  \alpha = 0.6 \pm 0.1.
\end{math}
These values could not be easily determined quantitatively, yet we successfully used qualitative observations to measure the filter's performance.

\end{document}
